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Abstract. Quite a few future high precision space missions for testing Special and General 
Relativity will use optical resonators which are used for laser frequency stabilization. These 
devices are used for carrying out tests of the isotropy of light (Michelson-Morley experiment) 
and of the universality of the gravitational redshift. As the resonator frequency not only 
depends on the speed of light but also on the resonator length, the quality of these 
measurements is very sensitive to elastic deformations of the optical resonator itself. As 
a consequence, a detailed knowledge about the deformations of the cavity is necessary. 
Therefore in this article we investigate the modeling of optical resonators in a space 
environment. Usually for simulation issues the Finite Element Method (FEM) is applied in 
order to investigate the influence of disturbances on the resonator measurements. However, for 
a careful control of the numerical quality of FEM simulations a comparison with an analytical 
solution of a simplified resonator model is beneficial. In this article we present an analytical 
solution for the problem of an elastic, isotropic, homogeneous free-flying cylinder in space 
under the influence of a tidal gravitational force. The solution is gained by solving the linear 
equations of elasticity for special boundary conditions. The applicability of using FEM codes 
for these simulations shall be verified through the comparison of the analytical solution with 
the results gained within the FEM code. 
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1. Motivation 

Special (SR) and General Relativity (GR) are two of the most important theories and 
theoretical frames of modern physics. They are the basis for the understanding of 
space and time and thus for the underlying physical structure of any other theory. 
The interest in testing the fundamentals of SR and GR has grown enormously over 
the last years as all presently discussed approaches to quantum gravity predict tiny 
violations of SR and GR. 

The technological improvements of the last decades have provided scientists 
with high precision measurement equipment such as optical resonators. In optical 
resonators laser locking is used to define stable optical frequencies. The resonance 
frequency of the locked lasers is given by v — mc/L where c is the speed of light, L 
the resonator length L, and m the mode number. 

Optical resonators have been used recently to test one of the pillars of Special 
Relativity, namely the isotropy of the speed of light [3 |H] as well as of the universality 
of the gravitational redshift 0. In doing so, two laser beams are locked to two 
orthogonally oriented resonators. An anisotropic speed of light would lead to a beat 
of the frequencies during a rotation of this setup. Due to the importance of this type 
of experiments one looks for ways to improve that. One option for this is to carry out 
these experiments in space, as planned in the OPTIS mission JH| or with SUMO on 
the ISS [IS|. 

Although many of the disturbances acting on a resonator can be minimized by 
means of an appropriate satellite control system, some intrinsic disturbances cannot 
be eliminated as a matter of principle and distort the resonator shape leading to a 
systematic frequency shift. In particular the tidal gravitational force \ which acts 
through every extended body cannot be eliminated by choosing an appropriate frame 
and, thus, will induce distortions on the resonator. 

We give a rough estimate of the expected effect of the tidal gravitational force on 
a freely moving cube of length L. If the position of the cube is at a distance R from 
the center of the Earth, then the difference of the Earth's acceleration on the top and 
bottom of the cube is Aa = (d 2 U/dr 2 )L, where U is the Earth's Newtonian potential 
U = GMq/R. For an orbit with R = 10000 km and a typical resonator length of 
L = 5 cm we have Aa w GAf®/ R 3 L w 2 ■ 10~ 8 m/s 2 . In a rough estimate we assume 
this Aa to act on the top surface of the cube. Now Hook's simple law of elasticity 

A L y ' 

gives the change of the length AL of the cube due to a force F acting on the area A. 
In our case taking F = mAa = pL 3 Aa we get 

^ = §Aa«10-- (2) 

assuming an elasticity modulus of E = 90 GPa and a density of 2350 kg/m 3 which is 
typical for Zerodur. 

In the OPTIS mission, for example, the science goal for the measurement of the 
isotropy of the speed of light is better than Ac/c = 10~ 18 This can only be 

achieved if the resonator has a length stability of AL/L = 10~ 18 |T^]. As one can see 
from our estimates, the tidal gravitational force will lead to systematic deformations 

X In space and engineering sciences the tidal gravitational force is often referred to as 'gravity 
gradient'. 
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which are one order of magnitude larger than the expected accuracy. Therefore the 
effect has to be investigated carefully by including the tidal gravitational force into 
the equations of elasticity, calculating the resulting resonator shape, and then subtract 
the effect. 

Although the linear theory of elasticity has a long history, explicit solutions for 
special problems are purely spread. In textbooks only examples for simple bodies 
in homogeneous gravitational fields or for thermal expansions can be found (e.g. 
|14| . [15| . |18|). However, most of the solutions employ an ansatz which already 
includes knowledge about the expected solution. To the understanding of the authors, 
no publications are available dealing with a body under the influence of a tidal 
gravitational force so far. The reason for this is probably, that this situation applies 
only to bodies freely flying in space - a situation which was outside the scope of 
application in elasticity theory so far. 

In the present paper we first derive an analytical solution in terms of a series 
expansion. This result is then confirmed using numerical methods. These calculations 
are usually done with help of Finite Element Method (FEM) codes. For most 
engineering purposes FEM codes are fine. However FEM solutions are only numerical 
approximations whose accuracy depends highly on the number and shape of the 
elements that have been chosen to mesh the model. In order to confirm the analytical 
model and to test the numerical calculation, we compare the analytical with the 
numerical solution. For this comparison we choose a cylinder as most simple geometry 
of a body adapted to the symmetry of the problem. 

Having thus checked the principal applicability of the FEM methods to these kinds 
of physical situations, this method safely can be used for calculating the deformations 
of arbitrarily shaped bodies or for the design of devices insensitive to unwanted 
influences, or for the elimination of the systematics of the measurements in order 
to ensure the success of highly sensitive experiments. 

2. Basic Equations 

2.1. Generalities 

The problem of an optical resonator flying on a geodetic Earth orbit can be simplified 
by treating the problem in a body fixed coordinate system. We also consider, for 
simplicity, the body to be a homogeneous and isotropic cylinder. The only force 
present is a volume force due to the tidal gravitational force which will be modeled as 
gradient of a spherically symmetric Earth acceleration field. 

In order to calculate the elastic deformations of the cylinder the equations of 
elasticity have to be solved including the influence of the tidal gravitational force. The 
boundary conditions for the solution are given through the condition of weightlessness 
in space. 

As a short introduction, some basic equations of the linear theory of elasticity are 
given |15M14I ITHlll 1| . All equations refer to homogeneous isotropic bodies. Within this 
paper we do not use the notation within the formalism of the Riemannian geometry 
(e.g. 120]) but the notation used in [T5j . 

In elasticity the general relation between the stress tensor o~ij and the strain tensor 
Eij is given by Hooke's law 

o~ij = Cijki£ij k,l — 1, 2, 3 (3) 
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where Cijki is the elasticity tensor related to the material under consideration. For 
homogeneous isotropic materials the elasticity tensor can be written as 

Cijki = A(5y Ski + n(SikSji + SaSjk) , (4) 

where A and fi are the Lame constants and b mn is the Kronecker symbol. Thus Hook's 
law for homogeneous isotropic materials is 

(Jij = XSijEkk + 2fJ,eij . (5) 

The strain tensor e has to fulfill the so-called compatibility condition 

d 2 e l i | d 2 e jk d 2 Eji d 2 e ik _ Q 



drjdrk dridri dridr k drjdri 

where r n are the components of the position vector. The relations between strain and 
the displacement & are 

The equilibrium equation of elasticity describes the equilibrium state of a 
homogeneous isotropic body when a volume force K is acting 

^a ij + K i = 0. (8) 

Applying the relations between stress and displacements the equilibrium equation 
takes the form 15 

fi + ri^b + ^fi+K-O, (9) 
where £ is the displacement vector. This equation can also be written as 

A| + ■ |) + -K = (10) 

where v is the Poisson number which lies between and 0.5 for homogeneous isotropic 
bodies. 

For vanishing volume forces K = Eq. 10 becomes the homogeneous equilibrium 
equation 

The boundary conditions for the solution of the equilibrium equation are either 
given by the forces pi acting on the body surfaces 

o-ijUj = pi (12) 

or by initial displacements £io of the surfaces 

6(0) =&> (13) 

where rij are the normal vectors on the surfaces. 

The general solution of Eq. I|10l) is a superposition of a homogeneous and a 
particular solution 

(14) 
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2. 2. The symmetries of our problem 

Since we have an axial symmetric problem, we use cylindrical coordinates r, ip, z is 
useful. All displacements and derivatives with respect to ip vanish and the equilibrium 
equation of elasticity takes the form (see e.g. |15p 



2 1 — 2v dr \ dr r dz 
Id fd£ r Cr d£A , . 

o = A6 + -^-^he + - + 7e • (15) 



1 — 2v dz \ dr r dz 
The Laplace operator acting on a scalar takes the form 
d 2 Id d 2 

dr 2 r dr dz 2 ^ ^ 

Note that the Laplace operators acting on a vector held £ takes the form (see |12) 1 

d 2 ir , ld£ r , 9 2 £ r £ r 



A£ r 



<9r 2 r 9r <9z 2 

9 2 6 , , a 2 ^ 



The relations between stresses, strains and displacements are 



+ e zz ) + 2fj,e zz = A 
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0~ r <p = 
G<pz = 



<t„ = 2 M e rz - /i ^— + — j . (18) 

Beside the axial symmetry we also have the following symmetries for reflection at 
the z = plane: £z(r, —z) = — £ z (r, z) and £ r (r, — z) = £r(f*, 2). 



3. The Problem 



In order to solve the problem of a free-flying isotropic homogeneous cylinder in space 
the equilibrium equation of elasticity (|10|l has to be solved. The cylinder has radius 
R and height 2L. The body coordinates are (r, ip, z) with the origin being at the 
center-of-mass of the cylinder. The z-axis coincides with the symmetry axis of the 
cylinder, see Fig. ^ 

The only force present is the volume force K which is due to the Earth's 
gravitational potential U, 

A|* + -J— V(V • - -pWU = . (19) 
1 — 2^ u 
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Figure 1. Simplified model of an optical resonator on a geodetic Earth orbit 



For a spherical Earth potential, U(r) — GM^/r, where GM§ is the gravitational 
constant times the mass of Earth, the potential acting at an arbitrary point P inside 
the cylinder can be calculated via Taylor expansion 

U[fM+rl = u(ni)+ ^ n+ iq^i nr , m 

CM 

= U(f M ) + WU(r M )r + ^(r 2 - 2z 2 ) (21) 

where ru is the vector from the center-of-mass of the Earth to the center-of-mass 
of the cylinder and r is the vector from the cylinder center-of-mass to point P. This 
Taylor expansion around the center-of-mass of the cylinder to second order gives the 
axis-symmetric potential in cylindrical coordinates. Note that the linear term of the 
Taylor expansion vanishes as this equation is valid in the freely falling reference frame 
of the cylinder. 

Since we consider a freely flying cylinder in an orbit around the Earth, no external 
forces are present and, thus, the forces p at the cylinder surfaces are zero which gives 
us the boundary conditions 

aijUj =pi = . (22) 

The normal vector n = {n r n v n z ) T (T means the transposed vector) reduces in 
the axis-symmetric case to n = (n r 0n z ) T , as the <p component is zero. 

Thus the boundary conditions at the top surface of the cylinder, i.e. for 
z = L,n= (001) T , are 

p r (r, z = L) = = a rz (r, z = L) 

Pz {r,z = L) = = <j zz (r 7 z = L) . (23) 
At the bottom surface of the cylinder, i.e. for z = — L, ft = (0 — 1) T , we have 
p r (r, z = -L) = = -<7 rz (r, z = —L) 

Pz (r,z = -L) = 0=-a zz (r,z = -L) . (24) 
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Note that these boundary conditions are valid for arbitrary r 6 [0, R). For r = R 
they are not valid as the normal vector is not uniquely defined at the cylinder edges 
(r,z) = (R,±L). 

For the superficies cylinder surface, i.e. r = R, n — (10 0) T , the boundary 
conditions are 

p r (r — R, z) = = a rr {r = R, z) 

Pz (r = R,z) = = cr zr (r = R,z) (25) 
for all z e [0,±L). 

4. The Solution 

4-1- Particular Solution 

In order to find a particular solution of the problem one can assume that the solution 
of the equilibrium equation can be written as gradient of a scalar tp |18| 

|p = . (26) 
Inserting this approach into Eq. I|19|) yields 

V f Ai/j + — - — Aip) = -VU 
\ 1 — 2// / /i 

K^ a *) =>■ (27) 

Calculation of the volume integral of the divergence gives 

A^^_p (C / + C / ) (28) 

where Uq is the integration constant. Inserting the potential U calculated in Eq. (|20|l 
and using (|16fl we obtain 

1 - 2v GAf ffi / r 4 1 4 2,2 , ,\ , nns 

2(1 — i/j^i 2r° M \ 16 6 / 
where Uq ~ Ac + 2d. Thus, the displacement vector becomes 

? - -5^^ ((-t + 2CT + *) * + (r 3 + M = + *) *) -< 30 > 

Since the displacement at the center of mass of the cylinder should vanish, £ r (0, 0) = 
£z(0, 0) = 0, both, c\ and d\, have to be zero. The unknown constants c and d are 
arbitrary. 

The r and z components of the displacement are (with p = A v" ) 
C P =-7T7^ rr/>T^ --T + 2cr = 77; T7 -- r + 2cr (31) 



2(1-1/)/^ 2< V 4 y A(l-iz) 

1 - 2v GMm (2 , „ l \ v (2 

— — p — ~z +2dz) = — -7 -i 

2(1 -v)^ 2r\t \3 ) A(l-zz) 1 \3 



where we substituted 7 := -GM/(2r 3 M )p. 
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Using Eqs. (|18|l we obtain the stress components out of the displacements 

,4(1-1/) \-v \-v 

^ CV V „2 , o^2 



1-1/ 1-1/ 



Homogeneous Solution 



r z + 2(V + d) 



(32) 



Jf.. 2.1. Derivation of the boundary conditions A displacement vector which satisfies 
the homogeneous equation of elasticity 

Af + y^;V(V • f ) = (33) 

has to fulfill the biharmonic equation (see, e.g. JH] and 0]) 

V 2 V 2£fa = AA £h = . (34) 

Note that the Laplace operator applied to a vector field has the form given in Equation 

03. 

The boundary conditions the homogeneous solution part has to fulfill can be 
derived from the boundary conditions (|22H of the complete solution. The boundary 
conditions at the cylinder top surface, i.e. z = L, normal vector n = (0 1) T , are 

p r (r, z = L) = = cr rz (r, z = L) = cr rz (r, z = L) + a^! z (r, z = L) 

p z (r, z = L) = = a zz {r, z = L) = o* z {r, z = L) + a p zz {r, z = L) . (35) 

With help of the expressions of the stress components derived from the particular 
solution (1M2I) one obtains 

(r, z = L) = — <(r, z = L) = - 7 -^r 2 + 2(L 2 + d)) 

\ 1 — v 1 — v J 

o* z (r, z = L)= - a p rz (r, z = L) = . (36) 

The boundary conditions at the cylinder bottom surface, i.e. z = —L, normal 
vector, n = (0 0— 1) T 

p r (r, z = —L) = = -cr rz (r, z = -L) 

= — o- rz (r, z = —L) — o-% z (r, z = —L) 
p z (r, z = —L) = = -o- zz (r, z = —L) 

= - a h zz (r, z = —L) — < (r, z = -L) (37) 
give the same boundary conditions as for the top surface, as z occurs only as z 2 term 
<j h zz (r,z = -L) = -o» z (r,z = -L) 

-r 2 + 2{L 2 + d) 



1-1/ 1 

cr^(r, z = —L) = - o-V z {r, z = -L) = . (38) 
For the superficies surface, i.e. r = R, n = (1 0) T 
p r (r = R, z) = = oy r (r = R, z) = cr^ r (r = R,z) + a^. r {r = R, z) 
p z (r =R,z) = = a zr (r =R,z) = a\ r (r =R,z) + a* r {r = R, z) (39) 
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one gets 

o-^ r (r = R,z) = - crj? r (r = R, z) 

/(2i/ — 3) ~ 2{c + vd) 2u 2 

\4(1 -i/) l-i/ 1-1/ 

o* (r = R,z) = - < (r = J2, z) = . (40) 

4-2.2. General Ansatz for the Homogeneous Solution Love ^3 showed, that the 
deformations in an elastic axis-symmetric body can be expressed in terms of the so- 
called Love function \ 

th. 1 + v d 2 x 



E drdz 



where £ r and £ z are the displacement components. 
The corresponding stress components are 
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The Love function x necessarily fulfills the biharmonic equation 

V 2 V 2 x = . (43) 

The major obstacle is to find an adequate approach for the Love function fulfilling 
all boundary conditions. A general approach is the separation of variables, that means 
a factorization of the kind xi r , z) = F\(r) ■ F 2 (z). For axial symmetric problems 
the Bessel functions (see |Appendix A| ) are the natural choice to represent the r 
dependence. Owing to the additional reflection symmetry and antisymmetry of £ r 
and £z , the z dependence can be represented by the trigonometric functions cosine and 
sine. Bessel functions as well as trigonometric functions form a complete orthogonal 
set of functions so that each function can be represented as a series of these sets. 

For the solution of the homogeneous equation of elasticity a so-called Papkovich- 
Neuber approach for the displacement field can be made J^l- O ne writes £ = 

V$ + 4(1 - v)B — V \B ■ f + B ), where $ is a scalar, and the vector B fulfills 



the biharmonic equation. For B and Bq the approach of separation of variables is 
used where the r dependency is represented by Bessel functions and the z dependency 
is represented by trigonometric functions. By modifying this approach one can make 
an ansatz for the Love function suggested in |2J 



\ - / , sinh(Ajz) coshpvjz)^ J (Xj 



X = B z A + > I Aj . + Bjz- 



\ 3 sinh(AjL) J sinh(AjL) J X) 



OO 



n i , \ ■ , n h{k n r) h(k n r) \ sin(fc»z) 



n=l 
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Herein J are the Bessel functions of first kind and order zero and I\ are the modified 
Bessel functions of order one (see |Appendix A| |. Furthermore Q = XjR are the zeros 
of the Bessel functions of order one, Ji(Cj) = 0, and k n = ^ where n is an integer 
number. 

If we insert the Love function approach into l|42(l we obtain for the stress 
components 

o* r = 6uB Q + {Av - 2) A, (45) 

E f , , cosh(A 1 z) „ /. „ cosh(A 1 z) , sinMA.z) \ \ ^ , 
_ sim^ + *> ( (1 + 2 ^imhTA7iy + ^imhfcl ) ) M ^ 



cosh(Ajz) /cosh(Ajz) sinh(Ajz)\\ Ji(Xjr) 

h(XjL) + j ^sinh(Aj-Z 

I (k n r) 1 h{k n r) 



El CUSn^AjZj / COSWyAjZ) Hi 

_ x \ srnh{Ayi) + J \sinh(AjL) + jZ iinh(Aji) 7 ; A,r 

OO / 

Sl CA UM) k n rh(k n R) 

+Al C 7 " "7-7-77T + - 2 " ! 77T7-77T ' ( '° MA '" : 1 

o*, = (6 - 6i/)-B + (8 - 4^)£> 

cosh(Ajz) /. . cosh(Aj-z) sinh(Ajz) 



+ 1 (<^.m + °- (« 2 - 2 ^ + ■ 

The appearance of the coefficients A,- and Bj in the Dini and Bessel-Fourier series 
(<^o + Y^=\djJo{Xjr) and Y^nLi c j^i(^i r )j see |Appendix A| | allows us to represent 
arbitrary boundary conditions for a zz and a rz at the cylinder top and bottom surfaces 
(z = ±L). Similarly, the appearance of the coefficients C n and D n in the Fourier series 
(ao/2 + a »i cos(/c„z) + XmLi sin(fc n z), see I Appendix B| ) allows us to describe 

arbitrary boundary conditions for o~ rr and o~ rz at the superficies surface of the cylinder. 
By using the boundary conditions we can now determine the unknown coefficients. 

4-2.3. Determination of the Coefficients Aj and C n For the boundary conditions 
= cr^(r, z = ±L) (Eqs. Qg) and we get 

OO 

= ^ (Aj-Aj + A (2f + Aj-icoth(AjL))) Ji(Aj-r) . (46) 
i=i 

Multiplication with kJi(Xj>r)rdr and subsequent integration over the interval to R 
yields with help of (|A.14(I 

A/A,,- = -Bj {2u + XjL coth(Aj-L)) j = 1 . . . oo . (47) 
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For the boundary conditions = of z (r — R, z) (Eq. I39fl we get 



sin(fc„z) . (48) 



Multiplication with sin(k n i z)dz and subsequent integration over the interval — L to L 
yields with help of (|B.2|) 



C n k n = -D n ( (2 - 2v) + fe»-R ^7^ J n = l...oo. (49) 
V Ii{k n R)J 

4-2.4- Determination of the Coefficients Bj and D n A simpler expression for the 
boundary conditions can be obtained by inserting <{4T|) and ^4*51 into (@5j), 

a* r (r = R,z) = -aP r (r = R,z) 

= 6vB + {Av - 2)D 



E'° ( , , %% cosh(A.,z) sinh(A 1 z)\ 



3 = 1 

oc 



°%z(r, z = ±L)= - a% x (r, z = ±L) 

= (6 - 6v)B + (8 - 4v)D 



V B 3 ( coth(A^) + ) JoM 



3 = 1 

OC 



+ Y(-l)»D n ( fc„r#4 + (2- k n R Io [ knR l) J » (fc " r Q (51) 



71=1 

Again Ji(XjR) = and cos(fc n L) = cos(n7r) = (— 1)" were used. 

Expansion of the Particular Boundary Conditions in Dini and Fourier Series As the 
Bessel functions form a complete orthogonal set of functions, each function can be 
represented by a Bessel series as given in (|A.6|) . With the additional condition that 
Q = XjR are the zeros of J\ and J| Jo(Ajr)| r= A = — Xj Ji(Xjr)\ r —u = —XjJi(XjR) = 
—XjRJ\(Cj) = 0, the particular boundary condition part o% z (r, z = ±Z) can be 
represented by a Dini series expansion as given in Equation IjA.llJl 

00 

(r, z = ±L) = a\ z (r) = d + ]T d 3 J (Xjr) (52) 

3=1 

with 

d " = ~EP? \ a zz( r ) rdr = Co - -^7- (53) 



i? 2 . 

where c = 7 (yzt; + ^(£ 2 + an< ^ c 2 = ^T~^ were substituted and 



4 " 

-7- j = 1 . . . 00 . (55) 



A3Jb(C,-) 1-" 
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Note that, following l|A.16p . 1 rJo(Xjr)rdr — if Q — XjR are the zeros of J±. 
Therefore the particular boundary condition cr^ z (r, z = ±L) can be written as 

cr L( r ) = co - A E , 2 r 77 , c 2 J (Xjr) . (56) 

Similarly, each function can be represented by a Fourier series, Eq. (|B.3(1 



o-P r (r = i?, z) = a? r (z) = y + a n cos(fc„z) + ^ 6 n sin(fc„z) (57) 

n— 1 n— 1 

and the Fourier coefficients can be determined from (|B.4I) . This gives a new 
representation of the particular boundary condition c^ r (r = R, z) 

OO T 9 „ i r 9 * f i \n 

a o L c 2 AL C2 ( — 1) 

°fr( Z ) = y + 2^ °™ C°S(fc„z) = C\ -\ 1 ^ ^2 COS(fc„z) , (58) 

n— 1 n— 1 

where 

/(2i/-3)„j 2(c + i/d)N , A 2v 

61=7 171 ^ + 1 and C2=7 1 59 

\4(1 — v) \ — v J 1 — 1/ 

with the arbitrary constants c and <i from the particular solution part. 

Inserting the Particular Boundary Conditions After having found representations 
for the particular boundary conditions in terms of cos and Jo we can determine the 
remaining unknown coefficients from Eqs. (|50|l and H51(l . In doing so we first simplify 
these equations with help of ljB.5j) . (|B.6(I and (|B.T|) 

- o-V r (r =R,z) = o*.(r = R, z) (60) 

L 2 c AL 2 c °° (— l) n 
-ex - —1 =-- V cos(fc„z) = 6vB + (4w - 2)D 

n— 1 

+ E B i E iTj^TW cos(fcnZ) Jo(Aji?) 

j'=l n=l » n j' 

-a z P 2 (r,z = ±L) = a*(r,z = ±L) (61) 

-co + ^— + J] TTT77T 62 J °(V) = ( 6 " 6 ^ B ° + ( 8 - 4 ^ D ° 
4 JT^ A j J o(0) 



00 / \ r \ 

gM cotli(w+ i^) j ° (v) 

Jo(Xjr) 

- i?(fc 2 + A 2 ) 2 J (A,i?) 



E(-l)»A>*n ( E 

71 = 1 



Multiplication of (|60|l with cos(k n /z) and subsequent integration over the interval 
from -I to L gives with help of | |B~2*) | 

J— 1 J 
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This results in an equation for the coefficients D n in dependence of B n 



h{k a Kf J k„R 



\m — 1 



In the same way multiplication of 161|) with Jo(Aj/r)r and subsequent integration 
over the interval from to R gives with help of ljA.151) 

^c 2 = MX 3 R)B 3 (coth(A^) + . y 



Af uw 7 J \ K J ' sinh^Aj-L) 



+ B- 1 )"^(hmta^J ' (64) 

what gives an equation for the coefficients Bj in dependence of D n 
^=(j (A^)(coth(A,L) + -^)) \ 

*{Tf-^-v' D ^{Tmfhw)) ■ (65) 

4-2.5. Determination of B and D a After the determination of the coefficients in 
the infinite series (|44f> we determine the remaining unknowns Bq and Do. From the 
boundary condition a^ r (r — R, z) — —af r (r — R, z), valid for all z € [0, ±L), one can 
deduce that in particular c\ r {r = R, z = 0) = —o-^ r (r = R, z — 0) must hold. Then 
Eq. (jSOJ) simplifies to 



7r^ * — ' n 

n—l 



(1 - XjL coth(XjL)) 

h{k n R) 2 \ 2-2v 



^ Bj I suffi^ ) Jo(A ^ } 



_ x n \ n lip 5 7 



(66) 



From E n =i = - 12 we § et 

-ci =6isB Q + (4 : v~2)D Q + n, (67) 
where we abbreviated 

'(1 - A i icoth(A JJ L)) 



i=i 
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Similarly, from the boundary condition o~^ z (r, z = ±L) = — cj z (r, z = ±L) valid 
for all r E [0, ±i£), one can deduce that in particular o\ z (r = 0, z = L) = — a\ z (r = 
0, z = L) must be valid. Then Eq. I|51|) simplifies to 

- co + ^-c 2 + V ^ 2 2 c 2 = (6 - 6^)B + (8 - 4z/)D + Z , (69) 

where we abbreviated 

°° / \ T 

Z=J2 B j coth(A 3 i) ' 



3 = 1 



sinh (XjL) 



1) D,^2 *nR h{knB) ) Jl(knR) 



(70) 



With joT^Tp = — "1" (compare Eq. l |B.12| l) we obtain 

- co = (6 - Qu)B + (8 - Av)D + Z . (71) 
Addition of Eqs. (EJ) and JJIJ g ive s 

- c - Ci = 6(£ + D ) + S (72) 
with S = 1Z + Z from which we can determine Bq and Do. Eq. H72|) yields 

B = -~(c + Cl +5) -D . (73) 
6 

and by inserting (f?3"|l into (|6Tjl we get 

Do = ^— ) {(l-v)c 1 -v(co+S) + K). (74) 

The unknowns c and d from the particular solution part can be chosen arbitrarily. 
Their influence on the homogeneous solution part is restricted to Bo and Do and is 
compensated in the complete solution. 

4-2.6. Summary: Coefficients of Homogeneous Solution For the sake of clearness 
the equations for the determination of the homogeneous solution part as derived in 
Eqs. ignj), O, CD> <E3 and JHH> are summarized: 

AjXj = - Bj (2v + XjLcoth(XjL)) (75) 

C n k n = - D n ((2 - 2u) + knE ^^l ) (76) 
V h{k n R) J 



7i(fc„i?) 2 7 fc„i? 
4c 2 (-l)« , ^ 4A m fc 2 



A../- 



Dj = Jo(XjR) coth(AjL) + - — f x 



sinh 2 (A_,L) 
2 °° / 4A 2 



Do = i (-co - c x - 5) - D (79) 
6 

D = 2{l l +v) ((l-v)c 1 -v(c +S)+K) (80) 



Deformation of a Cylinder under Tidal Gravitational Forces 15 
with 

{ Acv , , A f(2v-Z) o 2{c+vd)\ 2v , . 

<*=l\— v +*{L 2 + *= 7 (^# + -L_ij 7^(81) 

and S = 1Z + Z. 

Now we have determined all unknown coefficients in Equation l|44f) . Then, 
following Eq. (|4"T|l we can directly calculate the homogeneous displacement. By adding 
the homogeneous and the particular solution we get the total displacement of the 
cylinder. 

4-2.7. Convergence of the Homogeneous Solution The infinite system ()77l78fl can be 
approximately solved by reducing it to a finite system, that is, by expanding the sums 
only to n — N and j — J. Then we have a system of N + J equations. By increasing 
the values of N and J one can improve the accuracy of the solution and find their 
limits. 

This approach is justified, since we can prove that these infinite sums do converge, 
that means, that the infinite system (|77I78() possesses a unique bounded solution. In 
doing so we have to use the theory of regular infinite systems as has been done by ^2] 
and EH- 

In the following we will show first that the equation system (|77I78|) is a regular 
infinite system and second that it possesses a unique solution. For the sake of simplicity 
we introduce the abbreviations 

*lH + 5ife) (83) 

X n = (84) 



R 



A, 



Yj = — Bj^-Jo(XjR) (85) 
and rewrite the system (|77I78|) as 



1 ~ 4fc T 2 t _ _~ j_4a 2 



4A? „ „ ^ , „ 1 2c 2 



The system (jHEl can be written in a combined form 

oo 

Zj = ^ ejjzi + bj , i= 1,2, ...,oo, (87) 
i=i 

with Z21-1 = Xi and Z21 = Y\. Thus the vector z alternating contains the terms X n 
and Yj from Equation H86(l . Then the matrix e must fulfill 

e 2m -i,2i-i = (88) 

e2m-l,2i = c m,l (89) 

e2m,2;-i = dm,; (90) 
e 2m ,2/ = (91) 
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with 

1 4kl 1 4A| 

c ^'-p- n {kl + \]Y and 777:/.H • A;:^ ' (92) 

An infinite system of the form (j87(l is called regular, if in each equation of this 
system the sum of the norms of the coefficients is smaller than one (see |l()j ) 

oo 

X>i,*|<l (i = 1, 2, . . .) - (93) 
i=i 

Due to our substitutions we have 

oc oo 

£1^11 = 5^1^,11 (i = 2m-l) (94) 
i=i z=i 
00 00 

5Zl e Ml = 51^1 (* = 2m) (95) 

so that the condition for regularity reads 

00 

5Zl c wl<l (96) 
1=1 

00 

X)|dm,i|<l. (97) 
i=i 

These conditions are satisfied for the system Ij77l78(l or (I86|) because we can 



calculate using (|B.10(I and JB 

2 -''"" Pn^(^ + A|) 2 Pn U(W J K 

If 2-2v\ 4 
= — P, 



1 /2-2^-4 N 



1 



Pi V ^ 



2 



2 + 2f 

= 1 - *„ with = 

4A7 l L ( A; /. 



L = 7T L TI^XW = 777x7 coth(A^) 



+ Qj Aj V ' J ' sinh 2 (A,L) Aji 

= -Qi- A 

2 

= 1 - ^ with ^ = -A- . (98) 

Obviously, > as the cylinder length L and the radius R are positive and the 
zeros of the Bcssel function of first order XjR are positive, too. Using properties of 
the Bessel functions, it can be shown that also P n > 0. Therefore no singularities 
can occur. Furthermore '5 and $ are positive definite series and none of their terms 
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becomes one for all values of n and j, thus 1 — ^ n and 1 — $j are always smaller than 
zero. Therefore, (|S7|l is regular. 

The final conclusion that the regular infinite system I|77I78|I has a unique bounded 
solution can be inferred from a theorem |1(J| , that a regular infinite system <|87|) whose 
free terms fe, fulfill the condition 



with 



1 



E 



ei,i\ 



(99) 



does have a bounded solution \zi\ < K, where K is an arbitrary constant. 

From (PSJl its obvious that pi = 1 - |1 - ^| = and ^ = 1 - |1 - <F, 
Thus the condition \bi\ < Kpi can be divided into two conditions 



and 




< K\* n \ - A'i 



2 + 2^ 



k n Pn 



< A 



2v 



k 2 P 



< A 



Since K\ and are arbitrary constants they can be chosen to be 
4c 2 



K, = 



2 + 2v 



and 



K 2 = c 2 



(100) 



(101) 



(102) 



which fulfills the above condition. One can choose K to be K = K\ = ^+2v 



as the 

Poisson number v lies between and 0.5 and therefore K\ > K 2 . Thus we have proven 
that the condition for the existence of a bounded solution is fulfilled. 



5. Comparison between Analytical and FEM Solution 

Having proven that the infinite equation system (|77I78() possesses a unique bounded 
solution, we can now numerically calculate this solution with arbitrary order of 
accuracy by expanding the infinite series to sufficiently large numbers N = J. 
Although the infinite series converge very quickly, an expansion to higher orders N = J 
still gives an improvement of accuracy which can be seen in Table ^ 

Fig. 121 shows the resulting total displacement field from the analytical solution 
of Equation (|14|) . The infinite sums in the analytical solution were expanded to 
N = J = 1700. 

Now we can compare our analytical solution for our simplified problem with the 
result obtained with help of a finite element code in order to verify its applicability. 



Table 1. Comparison of the analytical solution at point r = I = R, z = 2 = L for different 
expansion orders N = J of the infinite sums in equation system 1771 . 



N = J £r 6 

100 -6.471 -lO" 15 1.637 • KT 13 

1000 -6.233 -10- 15 1.629 • 10~ 13 

1700 -6.218 -10- 15 1.628 • KT 13 
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4 4, A A A A k ^ 1^ 



\ \ \ ^ 1> V / r J l/ 

\ \ \ 1 \ J/ ' ]/ <l i/ l/ 



1 1 



Figure 2. Analytical solution: Displacement field of a cylinder under influence of a spherical 
tidal gravitational force field. The displacements are plotted over the body coordinates r and 
2. The cylinder boundaries are at z = ±L = ±2 and r = R = 1. 
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Figure 3. Finite element solution: Deformation of the cylinder under influence of a spherical 
tidal gravitational force field. The deformation is scaled by a factor of 6 • 10 13 . Right: 
Deformed cylinder shape and original finite element mesh. Left: the scale shows the z 
displacements. 



The finite element analysis was done with the commercial FEM code ANSYS. The 
cylinder model of radius R = 1 and length L — 2R was divided into approximately 
110 000 hexahedron elements. Hexahedron elements allow the creation of a structured 
finite element mesh which ensures a high relative accuracy of the finite element 
solution. Fig. El shows the deformation of the cylinder in the tidal gravitational force 
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Table 2. Displacements £ r and £ z in r and z direction; Comparison between analytical 
and FEM solution. The cylinder boundaries are r = R = 1 and z = ±L = ±2. In 
the analytical solution the infinite scries were expanded to N = J = 1700. The constant 
7 : = -^&p = - 3 - 98 6 °°j 4 1 l n 5 il° 14 2350 which corresponds to a 7000 km Earth orbit. 









Analytical 




Numerical 




r 


z 






tn 




fr 


c. 




1 


2 


-6.218 ■ 


10" ib 


1.628 ■ 10" 


■13 


-6.202 ■ 10" ib 


1.625 • 10- 


-ia 


1 





-3.491 • 


1Q -14 







-3.486 ■ 10" 14 


-1.66 ■ 10" 


-23 


0.49507\/2 


1 


-2.059 • 


io- 14 


1.151 ■ 10~ 


13 


-2.062 ■ 10" 14 


1.150- 10- 


-13 


0.11888^ 


1 


-5.208 ■ 


nr 15 


1.183- 10- 


13 


-5.236 ■ 10" 15 


1.182 - 10- 


-13 



field as result of the FEM calculation. 

Table contains a quantitative comparison between the displacements from the 
analytical and FEM solution for special points of the cylinder. Analytical as well as 
FEM solution are completely symmetric, i.e. the norms of the displacements of two 
opposing points of the cylinder are equal. Please note that we have chosen a very fine 
mesh for the finite element solution in order to get high accuracy. 

Obviously, the FEM solution and the analytical solution agree very well. The 
small differences are based on the nature of the FEM analysis. The elements in which 
the cylinder is divided cannot be chosen to be infinite small, they are 'finite' and thus 
the FEM solution is a kind of summation over all elements of finite size instead of 
an integration where the limit to infinitesimal small element size can be performed. 
Furthermore during the FEM analysis at least three points must be fixed in order 
to prevent the cylinder from rigid body rotations. In the current FEM analysis the 
center of mass as well as four of the next nodes were fixed in order to assure the perfect 
symmetry of the FEM solution. 

6. Summary and Outlook 

An analytical solution for the problem of an elastic, isotropic, homogeneous freely- 
flying cylinder in space under the influence of a tidal gravitational force has been 
worked out. The motivation for this work was the verification of the use of FEM 
codes for modeling the deformations of optical resonators in high precision space 
experiments. It could be shown that the analytical and the FEM solutions are 
equal within the limits set by the nature of the different solution approaches. The 
analytical solution itself gives a new application field for the theory of elasticity as so 
far (according to the knowledge of the authors) no examples with tidal gravitational 
forces have been calculated. 

A further application case of such analytical solutions is the estimation of the 
noise due to thermoelastic deformations in gravitational wave detectors. In particular 
with respect to the upcoming LISA (Laser Interferometer Space Antenna) mission 
|S] these investigations are of great interest. Some groups already dealt with this 
problem, see [51 116j. However, the solutions are given for infinite half-spaces only or 
some approximations have to be included in order to satisfy all boundary conditions. 
Therefore we are currently working on an exact analytical solution for the estimation 
of thermoelastic noise in gravitational wave detectors. 
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Appendix A. Bessel function relations 

For convenience of the reader we collect some equations from pp, (32] used 
within this article. 

The Bessel functions J v {x) of first kind and i^th order 

J_L (-l) v rx\ v + 2v 
vlL (v + v + 1) \2J 

V—O v ' 

are solutions of the Bessel differential equation 

x 2 y" + xy' + (x 2 - v 2 )y = . (A.2) 
Functions with the argument Ix fulfill the differential equation 

x 2 y" + xy' + (l 2 x 2 - n 2 )y = . (A.3) 
The modified Bessel functions I v {x) of first kind and vth. order are defined as 

00 1 +2 

W=r^ fa )=£ (|) " usR. (A.4) 

The orthogonality relations for the Bessel functions can be derived as 

R 2 ,, M2 rR 



y (^+l(U)) $ m n = J Ju{lumf)Ju{lun.r)rdr (A. 5) 

for v > — 1, v E M. provided that £„ m = l um R and £„„ = l vn R are the mth respectively 
the nth zero of J v , i.e. J v (l vm R) = and J v {l lJn K) — and it is shown that they form 
a complete set, so every arbitrary function f(r) can be represented by a Bessel-Fourier 
series 

oc 

/('') = X! C vrnJv (jv-m-j^) ( A - 6 ) 

for < r < R, v > — 1. The coefficients can be determined via 

Cvm = P2T7 ' 77 \12 / f( r ) J » f^mTi) rdr ■ ( A7 ) 

R [J V +l\C,vm)r Jo V R' 

Under the condition that l vm is related to the mth zero C,i/m of "^Jf^fvm"^) via 
Cum — lumR-i i- e - ■§pJv(lvm r )\r=R = 0, a second orthogonality relation for Bessel 
functions is 

— I 1 - 75— ) (J^(Cr/m)) S mn = / J v (lv m r)Jv(lv n r)rdr (A.8) 



2 V Z 



2 



for v > —\v > —1, v G 
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Because of this second orthogonality relation one can represent an arbitrary 
function as expansion of a so-called Dini series 

oo 

f(r) = d vm J u (Ui-^) (A.9) 

ni—1 

2 



Note that in the case of a Dini series for v = an additional term do has to be 
added in the series expansion |22| . For Dini series expansions and v = one has 

oo 

fir) = d +Y, d m J (l m ^) (A.ll) 



m— 1 

with 

2 



R 



; 77 m / /(r)Jo (/m^) rdr (A.12) 



and 

d ° = Jpj o f(r>dr (A.13) 

where £ TO = Z m i? is the mth zero of J\. 

Provided that l m — £ m /_R, where £ m is the mth zero of J\ and R is the maximum 
value of r one can derive the special relations 

R 2 2 f R R 2 2 

-^-(Jo(Cm)) ^mn = / Ji(l m r) Ji(l n r)rdr = — ( J 2 (Cm)) i5mn , ( A - 14 ) 

i? 2 2 

— (Jo (Cm)) ^mn = / Jo(l m r)J (l n r)rdr , (A.15) 
^ Jo 

= / Ml m r)rdr . (A. 16) 



Appendix B. Relations of Fourier and Dini Series 

The orthogonality relations for Fourier series are 
sin(n'z) sin(nz)dz = 



7T<W 


n 


7^0 





n 


= 




n 




2tt 


n 


= 



cos(n'2:) cos(nz)dz = 

sin(n'z) cos(nz)dz = (B.l) 
for integer n and n', n' ^ 0. 





n 


^0 





n 


= 




n 




2L 


n 


= 
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We also have 

sin(fc„/ z ) sm(k n z)dz = 

-L 
cL 

cos(k n i z) cos(k n z)dz = 

-L 
cL 

sm(kn'z) cos{k n z)dz = (B-2) 

' -L 

for k n — nir/L and k n < = n'ir/L, n, n' are integer numbers, n' =/= 0. 

As the trigonometric functions form a complete orthogonal set, each function f(z) 
can be represented by a Fourier series 

oo oo 

f(z) = — + ^2 a n cos(fc„z) + ^2 b n sin(k n z) . (B.3) 



The Fourier coefficients are determined via 

fi ,„_ , 1 rL 



I f / 717T \ ! f / n7T \ 

a n = — J /(z)cos^ — zj dz, bn ^i]J f (z) sin l^—z J dz . (B.4) 



In the following we collect some Fourier and Dini series expansions used within 
the article [§], |21|. for k n = nn/L and Ji(XjR) = 0, integer n,j: 

, sinh(Az) <T , „ r ,, cosh(A^) _,,„ 4AA; 2 , _. 

A2 ^Ai +(1 - ALcoth(AL)) ^h^ - T,yr L{kl+ \, r ^(k n z) { B.5) 



n=l 

h{kr) [2 UkR) \ I (kr) " 4Aj J (A 3 -r) 
' h{kR)J h{kR) ^R(k 2 + X 2 ) 2 Jo(XjR) ( ' 



h{kR) \k 



Io(kr) _ 2 " 2fc JqCV) 

/i (fc-R) fci? ^ R(k 2 + A?) Jo(A,-iJ) 1 ' ' 

From these equations one can derive the following sums, by setting z = L and r = R: 

> 77^ = T coth(AL) + 2- — - - — B.8 

^(fc2+A 2 ) 2 AV sinh 2 (AL) ALy 

y 4 ^ ff2 A ^o(fcfi) 2 \ , 2RI (kR) 

f^(k 2 + X 2 ) 2 V h{kR) 2 ) kh{kR) [ ' 

If /c — > in Equation (|B.9|) one gets 

In particular, 

°° 1 i? 2 

gX(V?)Af = -T- (R12) 
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